Name: Collaborators:
Due Wednesday 21
We want to solve the linear system given by \begin{equation} Ax = b, \end{equation} where $A$ is a non-singular $n\times n$ matrix.
In this homework you will find the $LU$ factorization of the Matrix $A$, such that \begin{equation} L\cdot U = A, \end{equation} where $L$ is lower triangular and $U$ is upper triangular.
If the LU factors are known, then solving the linear system can be performed in two triangular solves: \begin{equation} Lp = b, \end{equation} and \begin{equation} Ux = p. \end{equation}
Q1: Write a function that computes the LU factorization of a non-singular $n\times n$ matrix $A$. The input of your function will be
and the output will be a tuple $(L,U)$, such that $L\cdot U = A$, where
Moreover, your function will raise an error if:
You can use the built-in function det(), to check if you matrix is singular.
In [ ]:
function luFactorization(A)
# input A: a nxn matrix
# output L: lower triangular matrix
# U: upper triangular matrix
return (L,U)
end
In [ ]:
# write your triangular solver here
function backwardSubstitution(U,b)
# input: U upper triangular matrix
# b vector
# output: x = U\b
return x
end
Q1c: Write a function that implements a lower triangular solve. the input of your function will be :
your function will output $x$, the solution to $Lx = b$
In [29]:
# write your triangular solver here
function forwardSubstitution(L,b)
# input: L lower triangular matrix
# b vector
# output: x = L\b
return x
end
Out[29]:
In [ ]:
function linearSolve(A,b)
return x
end
As for Gaussian Elimination, we can encounter small pivots that may have a negative impact on the accuracy of the factorization. In order to avoid that we use partial pivoting (row pivoting) to reorder the unknowns.
Q1: Write a function that computes the pivoted LU factorization of the non-singular $n\times n$ matrix $A$. The input of your function will be
and the output will be a tuple $(L,U,P)$, such that $L\cdot U = P\cdot A$, where
Moreover, your function will raise an error if:
You can use the built-in function det(), to check if you matrix is singular.
In [ ]:
function rowSwap!(A, i,j)
# input: A matrix
# i,j the row indexes to be swapped
n = size(A)[1]
# checking that i and j are within the permitted range
(i > n || j > n ) && error("Index out of range in the rowSwap")
# if the index are not the same
if i != j
buffer = A[i,:]
A[i,:] = A[j,:]
A[j,:] = buffer
end
end
In [ ]:
function luFactorizationPivoted(A)
# write your function
return (L,U,P)
end
Q2b. Using your triangular solvers, write a function that solves the system $Ax = b$ using a partially pivoted LU factorization.
In [30]:
function linearSolvePivoted(A,b)
return x
end
Out[30]:
The complexity of LU factorization is $\mathcal{O}(n^3)$ for a dense $n\times n$ matrix.
Fortunately, most of the matrices arising from discretization of partial differential equations are not dense. In this case we will solve the 2D Helmholtz equation given by
\begin{equation} \triangle u + \omega^2m u = f \end{equation}plus absorbing boundary conditions.
The Helmholtz equation models waves propagating in an acoustic medium (like air or water). You will download the extra file Homework_7_extra.jl in the same folder as this ipybn, and you will run the snippet below.
We will first define the number of points
In [ ]:
# Pkg.add("PyPlot")
using PyPlot
include("Homework_7_extra.jl")
nx = 200
ny = 200
N = nx*ny
h = 1/(nx-1)
We will define the slowness $ m $ as a matrix full of ones, and we will call the function HelmholtzMatrix
In [2]:
m = ones(nx,nx)
# adding a reflector
m[round(Integer,nx/3):end, :] = 2;
H = HelmholtzMatrix(m,nx,ny, 10, h,300, 2,2*pi*nx/10 );
println("Average number of nonzeros per row is : ", nnz(H)/N)
As you can see, the matrix is extremely sparse, it only contains around 5 non zeros per row.
We will define a source, which in this case is a point source in the middle of the domain.
In [7]:
f = zeros(N,1);
f[nx*ny/2+ny/2] = 1/h^2;
Now we are ready to compute our solution. For that we will use the built-in LU factorization in Julia. (you may note that H is sparse, then by default julia will use UMFPACK to perform the factorization)
In [8]:
# perform the factorization
@time LU = lufact(H);
# solve the linear system
@time u = LU\f;
The solution is just a long vector, we need to reshape in order to plot it
In [5]:
U = reshape(u,nx,ny);
We can plot the solution using Gadfly, the solution is a complex number so we will plot the real and the imaginary part separately.
In [9]:
imshow(real(U))
Out[9]:
Q3a. You will benchmark the time it takes to compute the LU factorization and to perform the solve with respect to the size of the problem $N$. What is the complexity of the factorization and the triangular solves in this case? How is that different from what you have seen in class? How can you explain the difference?
In [ ]:
Hint: You may wnat to use the LU factorization algorithm you wrote (for a small $nx$) to obtain the LU factors. You can count the number of non-zeros of $L$ given by your algorithm and the number of nonzeros of $LU[:L]$, and compare them. You can use the function nnz(A) to obtain the number of nonzeros of a matrix.
In [ ]:
Q4a. Write a function that implements the fully pivoted LU factorization of the non-singular $n\times n$ matrix $A$. The input of your function will be
and the output will be a tuple $(L,U,P,Q)$, such that $L\cdot U = P\cdot A \cdot Q$, where
Moreover, your function will raise an error if:
You can use the built-in function det(), to check if you matrix is singular.
In [ ]:
function luFactorizationFullPivoted(A)
# checks
return (L,U,P,Q)
end